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Abstract. Analytical calculations based on finite-size spin-wave theory and Monte Carlo (MC) 
simulations are performed to investigate the validity of the well-known relation m(H, T) = 
M(H, T)Bd[M(H , T)J\fH/T] between the induced magnetization m of the magnetic particle and its intrin- 
sic magnetization M for the Ising and isotropic classical models \Bx){x) is the Langevin function, D is the 
number of spin components, A/" is the number of atoms in the particle] . It follows from general arguments 
and from our analytical results for the Heisenberg model at T < T c that this relation is not exact for any 
finite D and nonzero temperature. Nevertheless, corrections to this formula remain very small practically 
in the whole range T < T c if M S> 1, as confirmed by our Monte Carlo calculations. At T < T c /4 there 
is a good agreement between the MC and finite-size spin-wave calculations for the field dependence of m 
and M for the Heisenberg model with free boundary conditions. 

PACS. 75.50.Tt Fine particle systems - 75.10.Hk Classical spin models 



1 Introduction 

For magnetic particles of a finite size one can generally 
define two magnetizations, m and M, the relation between 
which is frequently written in the form 



m = MB D (Mx) , x = AfH/T, 



(1) 



where Bn{x) is the Langevin function [Bu(x) = B$(x) = 
coth x—l/x for the isotropic Heisenberg model and Bjj (x) = 
Bi(x) = tanhx for the Ising model] and AT is the number 
of magnetic atoms in the system. Here m is the magneti- 
zation induced by the magnetic field and microscopically 
defined as the thermodynamic average of the vector 



M 



i.e., 



rn (M). 



(3) 



For classical systems discussed throughout this paper, Sj 
can be considered, up to a factor, as spin vectors of unit 
length, | Si | = 1. The magnetization M in Eq. (Q) can be 
interpreted as the intrinsic magnetization of the particle 
which is defined through the correlation function of the 
magnetic moments, 



M = v/(M 2 ) 



(4) 



If the temperature is low, all spins in the particle are 
bound together by the exchange interaction and M be- 
haves as a rigid "giant spin", |M| =M=l,whicli shows a 
superparamagnetic behavior. If a magnetic field H is ap- 
plied, M exibits an average in the direction of H, which 
leads to a nonzero value of the induced magnetization 
m given, obviously, by Eq. ([l]). The question of princi- 
pal interest is, however, the field dependence of M at 
nonzero temperatures, which can be responsible for de- 
viations from the simple superparamagnetic behaviour of 
Eq.©. 

Early Monte Carlo (MC) simulations by Wildpaner 
for the classical Heisenberg model, where both magne- 
tizations were determined independently as functions of 
field at different temperatures, confirmed Eq. (Q) within 
numerical errors. However, from the theoretical point of 
view this relation with M = M(H,T) is unexpected. 

On the theoretical side, Eq. (|l]) was obtained in Ref. || 
for a classical model and in Ref. [[3) for a quantum model 
but without the field dependence of M. Earlier, Fisher and 
Privman Q considered the spin-wave contribution to Eq. 
(|l|) but, again, the field dependence of M was not studied 
explicitly. 

Experimentally, the field dependence of M and, in par- 
ticular, the nonsaturation of the magnetization in the re- 
gion x > 1 have been observed in nanoparticles by dif- 
ferent groups Usually this dependence is close to 
linear and is used to extract the value of M at zero field 
by extrapolation to H = 0. For the isotropic Heisenberg 
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model, the field dependence of M in the range x > 1 is due 
to suppression of the fluctuations of individual spins, i.e., 
of spin waves, and this dependence disappears for T — > 0. 
The dependence M{H) is much stronger and persists at 
zero temperatures if the spins in the particle are not per- 
fectly collinear due to surface effects 0. 

In our recent paper § (see also Ref.§) we have shown 
that this relation becomes exact for the exactly solvable 
model of the D-component classical vector "spins" in the 
limit D — > oo. Nevertheless, for more realistic models such 
as the classical Heisenberg model (D = 3) and the Ising 
model {D = 1), it is very difficult to believe that the super- 
paramagnetic relation holds for all temperatures. Clearly, 
if the number J\f of atoms in the particle is large and 
the temperature is below T c , then the argument of the 
Langevin function in Eq. (Q) becomes large already for so 
small fields that M does not essentially deviate from its 
zero- field value. Under these conditions Eq. (|I|) should be a 
good approximation. On the other hand, for smaller parti- 
cles and near or above T c , there should be deviations from 
the simple behavior, the study of which is the purpose of 
this work. 

The structure of the rest of this article is as follows. 
In Sec. H using the low-field expansion of m and general 
arguments we show that Eq. (|l|) is not exact for any finite 
value of D and nonzero temperatures. In particular, in the 
high-temperature limit there is another analytic form of 
Eq. (0) with Bp substituted by B^. In Sec. || we present an 
explicit calculation of both m(H,T) and M(H,T) at low 
temperatures with the help of a spin-wave theory which 
separates the global-rotation mode and the k ^ spin- 
wave modes. In Sec.^J we perform high-accuracy MC sim- 
ulations for the Ising and classical Heisenberg models in 
the box geometry to illustrate the superparamagnetic be- 
havior in a wide range of parameters. 



2 Basic Relations 

We use the classical spin- vector Hamiltonian 



n 



(5) 



where s is a £>-component vector (D = 1 for the Ising 
model and D = 3 for the Heisenberg model). For this 
Hamiltonian one can prove an identity relating correla- 
tions functions and susceptibilities 



M = m 1 



dm (D - l)r 



dx 



(0) 



where x is given by Eq. ([!]). On the right-hand side of 
Eq. (^) , the second and third terms are contributions from 
the longitudinal and D — 1 transverse susceptibilities, re- 
spectively. This relation can be used to extract the value 
of M(H,T) from measurements of the induced magneti- 
zation m and susceptibilities. Let us demonstrate how it 



works at low fields, where the expansion of m can be writ- 
ten as 

„2 „4 

-3 



m = —x ~ 



D D 2 (D + 2 
Applying Eq. (|^) one readily obtains 

M = < 



■X 



2aD 2 



(7) 



(8) 



At zero temperature the magnetic moment of the particle 
can be considered as a rigid spin, thus in Eq. (Rl) a = c = 1 
which results in M = 1, independently of the field. At 
T > one has a < 1 and c < a, so that M increases 
with the field. The coefficients a and c can be calculated 
analytically at low and high temperatures [see, e.g., Eq. 
(pOf)], Let us check now what happens if we try to find M 
from Eq. (m) under the same conditions. One can write 



M 4 



M 2 



D D 2 {D + 2) 



M ^ M Q + M 2 x 2 , (9) 



and find the coefficients Mq and Mi from the condition 
that m here coinsides with that of Eq. ^7§. The result is 
Eq. (||) with D 2 — * D(D + 2). This is clearly a wrong result 
for any finite value of D and nonzero temperature. Only 
in the limit T — > the coefficient M 2 vanishes and both 
approaches yield the same trivial result. Therefore, one 
cannot use Eq. ([!]) to take into account the field variation 
of M in the range where the argument of B is of order one 
or less. This formula can only be correct in the case of large 
particles for which the change of M in this field range is 
very small and M actually changes for much larger fields 
where we already have m = M. 

On the other hand, using these results one can find 
the correction to Eq. (|f|) at low fields. To this end, one can 
write 

m = MB(Mx) + S, (10) 

expand it for x <C I using Eq. (||) and equate the result to 
Eq. (0) to find S. The result is 



«J = - 



2(a 4 - c 4 ) x 3 
D + 2 TP 



<0, 



(11) 



that is, the Langevin function Bd in Eq. (|I|) should be 
replaced by some function F which goes below B^ at 
nonzero temperatures. 

In the high-temperature limit one can find an explicit 
form of the superparamagnetic relation which also differs 
from Eq. ([!]). Indeed, at high temperatures the exchange 
interaction can be neglected and one has to solve a one- 
spin problem, which yields 



M 2 = m 2 



1 



B'(0 + (£>-!) 



(12) 

where £ = H/T. Using this relation, one can plot m/M 
vs xm = Mx and thus obtain the scaling function F(x) 
which replaces Bp(x) in Eq. ([I]). For large particles, JV ^> 
1, in the relevant region xm ~ 1 one has ( < 1 and the 
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second of Eqs.@, with the use of B'{£) ^ - VA 

simplifies to 



M = m + Dm /a 



(13) 



On the other hand, this relation holds in the large- D limit 
for all temperatures, particle sizes, and types of boundary 
conditions, and it can be obtained from Eq. (|^) by drop- 
ping the term dm/dx and replacing D — 1 — > D. Solv- 
ing this equation for m yields the scaling function of the 
spherical model 



F(x)=B 00 (x) = 



2x/D 



l+y/l + (2x/D) 2 
in Eq. (|l|) which goes below Bjj{x) for any finite D. 



(14) 



3 Spin-wave theory for finite-size magnetic 
particles 

3.1 General 

At low temperatures all spins in the particle are strongly 
correlated and they form a "giant spin" M [see Eq. (g)] 
which behaves superparamagnetically. In addition, there 
are internal spin-wave excitations in the particle which 
are responsible at nonzero temperatures for the fact that 
M < 1 and for the field dependence of M. In our case 
of three-dimensional particles, d — 3, these excitations 
can be described perturbatively in small deviations of in- 
dividual spins Si from the direction of M. To this end, 
it is convenient to insert an additional integration over 
dM. = M dMdn in the partition function, 



-n/T 



Z = J M^dMdnJJdmS f M -^X) s ^ 

(15) 

and first integrate over the magnitude M of the central 
spin [this variable appears locally and it should not be 
confused with the intrinsic magnetization M defined by 
Eq. (Q)]. To do this, one should reexpress the vector argu- 
ment of the 5-function in the coordinate system specified 
by the direction of the central spin n, which yields 

x*f^£[si-n(n.Bi)]^16) 



Then after integration over M one obtains 



Z = dwZry 



where 



(17) 



(18) 



and 



n cff = -(n • H) ^(n ■ - i ^ J y" s * ' s o 



(D - l)Tln 



(19) 



In Eq. (|18|), the ^-function expresses the obvious condi- 
tion that the sum of all spins does not have a component 
perpendicular to the central spin M. This will lead to the 
absence of the zero Fourier component of the transverse 
fluctuations of spins in the particle. The corresponding 
global-rotation Goldstone mode (which is troublesome in 
the standard spin- wave theory for finite systems) has been 
transformed into the integration over the global variable n 
in Eq. ( |l7j ) in the present formalism. The condition men- 
tioned above was also used to transform the Zeeman term 
in Eq. (|l9|). This describes now the spins Sj in a field in 
the direction n and with the strength n • H. As we will 
see below, the last term in Eq. ( |l9| ) is nonessential in the 
leading approximation at low temperatures. 

To calculate Z n at low temperatures, one can expand 
7i c ff up to the bilinear terms in the transverse spin com- 
ponents 

IT = Si n(n • Si ) (20) 



using 



This yields 



n • s, ; = 



= Vi-n? = i-n?/2. 



(21) 



H cS Si E a - Mn • H 



n, 



At 



(D- 1)T/M + n-U+J2 Ju 



Jij, (22) 



where E = —(1/2)^ - Jij is the zero-field ground-state 
energy. For the lattice sites inside the particle or for the 
model with periodic boundary conditions one has z = 2d, 
where d is the spatial dimension; for the sites on the 
boundaries Zi < 2d. Now Z n in Eq. (Q) takes on the form 



Z n 



exp 



-ffo+A/'n-H 
T 



D-l 



D-l 



nn^ 



i a — 1 



x6 ^ exp ^--L J2 A v U i ■ > ( 23 ) 

which after working out the Gaussian integral over 77" 
yields 



Z n = exp 



-E +Afn-Il 
T 



D-l 



(2nT) 



JV-1 



detA' 



(D-l)/2 



(24) 



where the matrix A'^ is obtained from Ay of Eq. (|22j) by 
elimination of 11^ on one of the Af lattice sites using the 
condition J2i n* = 0. 
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Eq. ( g4j ) is a general result which is valid for a parti- 
cle of arbitrary shape and for different types of exchange 
interaction and boundary conditions. In the following 
subsection we will consider particles of cubic shape with 
the nearest-neighbour interactions and free and periodic 
boundary conditions (fbc and pbc). 



3.2 Free and periodic boundary conditions 

Let us express the matrix A$j through its eigenfunctions 
fki as follows 

^ = £/w4feAi> (25) 
k 



where fki satisfy 



only systems with real eigenfunctions. In this case, inte- 
grating over (Af — 1)(D — 1) modes (Af — 1 fc-modes mul- 
tiplied by D — 1 transverse spin components) for Z n one 
obtains Eq. (53) with 



det^,=n 

fe 



Af ' 



(32) 



where the prime on the product means that the mode with 
k = is omitted. 

All the results above are still general. Now we will con- 
sider cubic-shaped particles with free and periodic bound- 
ary conditions. In the fbc case the matrix Aij has the form 



A, 



A 8i 



(33) 



(26) where Aq is given by Eq. (J30|) and 



and form an orthonormal and complete basis 

^ fkifk'i = $kk> , 5Z fkifki' = ■ (27) 



as 



where 



53 Aijiii Uj = 53 A k u l ■ n fe , 



life = 53 /feillj 



(28) 



(29) 



Now one can make the observation that in the set of 
eigenfunctions fki there is one which is independent of i 
and which can be conveniently ascribed to k = 0, i.e., 
/o = l/Af. This follows since 

53 Aij = Aq = (D — l)T/Af + n • H (30) 



is independent of j, Aq being the zero-k eigenvalue. Now 
one can see that tfQ^IIi) m Eq. (23) excludes integra- 
tion over the zero mode IIo in the new representation. 
Fluctuations of the components life with k 7^ yield mul- 
tiplicative contribitions to the partition function, so that 
one is left with the integrals over 77^ . If the eigenfunctions 
fki are real, one obtains 



dllk exp 



Ak 
2T 



1/2 



br ■ (31) 



If fki are complex, one has to integrate independently 
over the real and imaginary components x^ and y^ of 
77^ = + iy^ which gives 2-kT jAk. Complex eigenfunc- 
tions arise, however, only in the case of periodic boundary 
conditions where, as we shall see, one has to take into ac- 
count only a half of the k modes, which effectively restores 
the result of Eq. (|3ll). So let us consider for the moment 



A 



(x) 



(34) 



etc., j x are Kronecker symbols, and i x ,j x = 1 5 



In this basis, the sum over ij in Eq. (E3j) can be rewritten § 



,N. 

If j x — 1 or j x + 1 run out of the particle, the corresponding 



( x \ 

A\j is a discrete Laplace operator for the coordinate x, 
and the eigenvalue problem factorizes. The eigenfunctions 
are standing waves and they can be written in the form 



-l or di xi j x + 



i should be omitted. One can see that 



/ki fi x ,kx X fi y .ky X fi z ,k z i 

where for a = x, y, z 



(35) 



i a ,k a — 



■ cos[(i Q - l/2)fc Q ], 



k a =im a /N, n a = 0, 1,...,N-1. (36) 



For the eigenvalue Ak one obtains 

A k = A + J k - Jo- 



(37) 



In the case of periodic boundary conditions, one should 
drop the terms and Sj xj N & n d identify i x = N+l with 
i x = 1 in Eq. (|34|) . The eigenfunctions can be conveniently 
taken in the form of plane waves e _lkri with the wave 
vectors quantized as k a — 2?:n a /N, the eigenvalue A^ 
having the same form as in the fbc case. That is, the pbc 
and fbc models differ only by the quantization of the wave 
vector 



2irn a /N, pbc 
nn a /N, fbc ' 



0,1,...,7V-1 (38) 



where a — x,y, z. This subtle difference is responsible for 
much stronger thermal fluctuations in the fbc model due 
to surface effects, as we shall see below. 
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3.3 The partition function 



with n 



and 



Collecting the formulae obtained above, one can write 
down the expession for Z n in the form 



Aexp ( n ■ x + TV 



D - 1 



/(Gn) 



(39) 



where A is given by the first line of Eq. (|56j) below, 

x = 7VH/T (40) 
is the reduced field, the function /(G n ) is defined by 

G n 

m — 

k 



G n Al, 



(41) 



with Ak = Jk/ Jo = (cos k x + cos k y + cos k z )/d, and 
1 T 



G n = — = 1 - a n , a n = —r— (n • x+D - 1) < 1. 

1 + a n TV J 

Note that the angular dependence of 2 n is more com- 
plicated than that for rigid spins because of the internal 
spin- wave modes described by the last term in Eq. (fj9|). 
These SW modes have a gap accounted for by the first two 
terms in the denominator of Eq. ( |32| ) or the dimensionless 
parameter a n in Eq. (f42|). One contribition to the gap is 
due to the finite size of the particle and the other is due to 
the magnetic field. The latter depends on the orientation 
of the particle's magnetic moment n with respect to the 
field. 

The function /(G n ) of Eq. @ can be written as 



/(G n ) = /(l) - f 1 du 



Pn(u) 



where f{\) is a constant and 



(43) 



(44) 



is the lattice Green function. Since at low temperatures 
the argument G in the expressions above is close to 1, it 
is convenient to write 



P N {G) = Pjv(1) 



-y- 

(1 



(1 - G)\y 



(l-GA k )(l-A k ) 



W N -N(l-G)f P (l-G), 



(45) 



where Wn = P/v(l)- Here, if the linear size JV is not large, 
one can replace G — * 1 in the argument of the function 
fp. For iV ^> 1 the situation becomes more complicated 
since the wave vectors k come closer to the origin and a 
singularity is formed. For the system with free boundary 
conditions, the sum is dominated by k <C 1, so that Ak — 
1 — k 2 /(2d) and fp(y) has the form 



fp(y) 



(2d)' 



1 



^ (y + n 2 )n 2 



(46) 



y = 2d(l~G)(N/ir) 2 



(47) 



For y <C 1 one can set y = which yields fp(y) — /p(0) = 
Cfbc — 1.90, whereas for i/> 1 one can replace summation 
by integration and calculate the integral analytically. For 
the model with periodic boundary conditions, there are 
different contributions from different corners of the Bril- 
louin zone in Eq. (|45|) , and one obtains a more cumbersome 
analogue of Eq. (f4q). In practice, it is easier to compute 
fp from its definition in Eq. (fi5|). For three-dimensional 
cubic particles the limiting cases are (1 — G <SC 1) 



P N (G) S W N 



c N N{\ 



G, y>l, 



(48) 



where for large AT the value of Wiv at 
integral W = 1.51639 according to f 



sroaches the Watson 



A 



N 



W N -W 

w 



0.90 

jy"' 

91n(1.177V) 
27rPFiV 



pbc 



fbc 



(49) 



(notice the positive sign for the fbc model). For the sim- 
ple cubic lattice cq = (2/-7r)(3/2) 3 / 2 . The numerically ob- 
tained results for cat can, for N 3> 1, be fitted as 



cn 



0.384 - 1.05/7V, pbc 
1.90 - 1.17/JV, fbc 



(50) 



(see Fig.|l|). The square- root term in Eq. flis] ) describes the 
spin-wave singularity in the infinite system. From Eqs. 
([l?]) and (^) it follows that the crossover to the bulk 
spin-wave behavior occurs for the values of the reduced 
field x > xg ~ NJ /T which is much larger than the 
value x ~ xy = 1 corresponding to the suppression of the 
global rotation of the particle's magnetic moment. The 
actual crossover fields, in notations of Ref. [Q, are given 

by 



H v = « H s 



Jo 



2dAf 2 / 3 ' 



(51) 



that is, they are widely separated from each other in our 
case T/(N Jo) <C 1. Thus the result for the function /(G n ) 
of Eq. ([i3|) , which with the help of Eq. (^) can be written 



as 



f(G n ) = /(l) - a n W N + N / daaf P (y a ) (52) 

Jo 

with y a = 2da(N / it) 2 , can be simplified in different field 
ranges. 

For H <§C Hs one can replace fp(y a ) by cat to obtain, 
in Eq. ©, 

Af^r^/ (G n ) £* tf£^±f(X) - t( n • x + D - 1) 



+a(n-x + L> - l) 2 , 



(53) 
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2,0- 
1 

1,6- 
1,4- 
1,2- 
1,0- 
0.8- 
0.6- 
0.4- 
0,2- 
0,0- 




1.90 - 1.17/W 




0.384 - 1.05/W 



xB{[l-t + 2a(D - l)]x} 
+a[2x- (D- l)(B + xB')] 
= [l-t + a(D-l)} 

xB{[l~t + a(D-l)]x} + 2ax. (59) 

Expanding the expression for m for x ^ 1 leads to Eq. (0) 
with the explicit values of the parameters 



a 2 = [1 - t 



a(D-l)} 2 
t + a{D-l)]\ 



2aD, 



(60) 



0,4 



1/JV 



0.5 



Fig. 1. Finite-size effect on cjv = /p(0) [see Eq. ([tH])] for cubic 
systems with free and periodic boundary conditions. 



Note that in the region x 3> 1, where a rigid magnetic 
moment would saturate, m continues to increase linearly 
as m = 1 — t + a(D — 1) + 2ax. This is due to the field 
dependence of the intrinsic magnetization M. The latter 
can be calculated from Eq. (||) which leads to 



where 



t 



M = l — t + a(2D - 1) + 2axB{x). 



(61) 



D - 1 W N T 
~2 JT' 



{D - l)c N (T_ 
4JV 2 I Jo 



(54) 



are small parameters, a <C t <C 1. Since ax 2 ~ (H/Hs) 2 <C 
1, one can expand the partition function Z of Eqs. ( Jl7| ) and 
(p9|) with respect to the last term of Eq. (p3), which yields 



This formula decribes a crossover from the quadratic field 
dependence of M at low field, x <C 1, to the linear depen- 
dence at x 1. 

Now we are in a position to calculate the correction 
to Eq. (fy) at low temperatures and x ~ 1. To this end, 



Z^A' 



one can write m in the form of Eq. (10), expand it with 
respect toa<l and equate to the expanded form of Eq. 
(0). This gives 



x exp { [1 - t + 2a(D - 1)] ux} [l + a{uxf] (55) 



6 = a[2x — (D + 2xB){B + xB')\ , 



(62) 



where 



A' = M - 1 (2nUT/Jo) iD - 1] ^- 1)/2 e- Bo/T 



xe (D-l)[A/7(l)/2-t] 5 



D-l 



(56) 



and Sd = 2n D I 2 j r(D j2) is the surface of the Z?-dimensional 
unit sphere. In fact, we have left the term proportional to 
aux in Eq. (^5|) not expanded for the sake of convenience. 
Integration in Eq. ( ps| ) results in 



Z = Z { [1 - t + 2a(D - l)]x} 



,d 2 Z Q 
dx 2 



(57) 



where Zq{ [1 — t + 2a(D — l)]x} is the partition function of 
the rigid magnetic moment with the magnitude reduced 
by the factor 1 - t + 2a(D - 1). 



3.4 The superparamagnetic relation 

Using 

1 d 2 Z n , (D-l)B 

= B + B =1 — ± '—, (58) 



which has a negative value. In particular, for x 1 one 
has S = -8ax 3 /[D 2 (D + 2)} [cf. Eq. @]. It can be shown 
that S — > in the large- D limit. Since a defined by Eq. 
( p4| ) contains N 2 in the denominator, 5 remains small even 
if T ~ Jo. This is an indication that Eq. (|l|) is a very 
good approximation for not extremely small systems in 
the whole range below T c . It can be shown that for ./V 3> 1 
crossover to the high-temperature form of Eq. (|l|) specified 
by the function B^ (x) of Eq. ( |l4| ) occurs in a close vicinity 
of T c . 

At higher fields H ~ H$ there is another crossover 
to the standard spin-wave-theory result for M. Here one 
has x ^> 1, thus the integral in Eq. ( |i~7j ) is dominated by 
n • x = x. Replacing n • x — > x in the last term of Eq. (52) 
one obtains 



Z = Zq[(1 — t)x]exp 



D - 1 



AT 4 / 3 / daaf P (y a ) 



a x = xT/(MJ ) = H/J 



Za dx 2 



which yields 
m = M = 1 



D — IT 
2 Jq 



W N - ^Nf P f-|- 
Jq \tis 



(63) 



(64) 



where B = B(x) is the Langevin function, for the induced 
magnetization m one obtains 

dlnZ , , „ 

m = — — ^ 1 - t + 2a(D - 1 
dx 



where the function fp(y) is defined by Eq. (|4q ) and H$ is 
defined by Eq. @. 

Let us now write down the explicit forms of the field 
dependence of the intrinsic magnetization M in the three 



H. Kachkachi and D. A.Garanin: Spin-wave theory for finite classical magnets and superparamagnetic relation 



7 



different field regions 



M =l — t 




1 L(t£\ 112 

_C ° Jo \ T ) 



H -C Hy 



H s < H < J . 



_ (65) 
Here i 1 is defined by Eq. (J54|) . In the second and third 
field ranges, the particle's magnetic moment is fully ori- 
ented by the field, thus m = M, the spin-wave gap in Eq. 
(32|) has the value H, as in the bulk, and the field depen- 
dence of both magnetizations follows that of the function 
Pjv(G) of Eqs. @ or @ with 1 — G = H/J [see Eq. 
(0)] • The region H <C Hy in Eq. ( |65| ) is less trivial. Here 
the gap in Eq. (^) is n • H and depends on the orienta- 
tion of the particle's magnetic moment which is not yet 
completely aligned by the field. Effectively one has in this 
region n • H ~£f 2 , which leads to a quadratic field depen- 
dence of M. In fact, such a dependence at smallest fields 
already follows from general principles, see Sec. ||, and is 
pertinent to the Ising model as well. 

To conclude this subsection, we introduce the orientation 
dependent "macroscopic" particle's magnetization M n ac- 
cording to 

d In Z n , „ 

M tl = — 66 
ax 

where Z n and x are defined by Eqs. ( |f~8| ) and (fiol), respec- 
tively. Using this definition, for the induced magnetization 
m = <91n Z/dx one can write 



which readily yields Eq. ( |6l| ) up to a field-independent 
term (a<i). 



3.5 Local magnetization 

The formalism developed above can be applied to study 
inhomogeneities in the particle's magnetization arising as 
a consequence of free boundaries. Since in zero field the 
standardly defined local induced magnetization = (sj) 
of a finite-size particle vanishes, one has to introduce local 
intrinsic magnetization 



f dn M n Z n 

m = 

J dnZ n 

M n can be interpreted as M of Ec 
modes integrated out. From Eq. (|3€ 



(67) 

m with the spin-wave 
one obtains 



1 



D-l T P N {G n ) 
2 Jq G n 



which for H <C Hs can be written as 



[I - t + 2a(n • x + D - 1)] n. 



(68) 



(69) 



The magnitude of the particle's magnetization, M n = 
|M n |, depends on its orientation due to spin- wave effects. 
It attains its maximal value I — t + 2a(x + D — 1) if the 
particle's magnetization is directed along the field H and 
its minimal value 1 — t + 2a(— x + D — I) in the thermody- 
namically unfavorable state with magnetization against 
the field. It should be stressed that in order to obtain 
the explicit result for the induced magnetization, Eq. , 
from Eq. (p^), one should know i? n , so its calculation in 
the main part of this section is unavoidable. On the other 
hand, for the intrinsic magnetization M it is sufficient to 
replace Z n =>• cxp(n • x) and use 



M 



J dnM n exp(n ■ x) 
J Jnexp(n ■ x) 



(70) 



(71) 



One can check the identity (1/Af) = M showing 

the self-consistency of the definition given above. Adding 
the expression within the brackets to the integrand of Eq. 
( |l5| ) and repeating all operations, one arrives at the final 
low-temperature result 



Mi 



D-1T 



dnZ„ 



k 



\fki\ 



(72) 



where and fki are eigenvalues and eigenfunctions of the 
linear problem, see Eq. (p6|). The latter contain the infor- 
mation about inhomogeneities in the system. For periodic 
boundary conditions, one has fki = e- lUr '/y/N, so that 
l/fcil = yl/iV and there are no inhomogeneities. Since the 
parameter t defined by Eq. (|54| ) is small, one can expand 
Eq. (|7^) to obtain, to the lowest order at low temperatures, 



Mi = I - 



(D — 1)T 



E 



a-; 



At 



(73) 



Here one can check again (1/AQ J2i % — M, where M = 
I — t, according to Eq. (61). For cubic particles with free 



boundary conditions, one has fki ~ \/2jN at the bound- 
ary according to Eq. ([36]), which is larger than the bulk- 
averaged value. The biggest effect of the surface is natu- 
rally attained at the corners of the cube where M, w 1— 8t. 



4 MC simulations and results 

The classical Monte Carlo (MC) method based on the 
Metropolis algorithm is now a standard technique (see, 
e.g., Ref. for details). The general idea is to simu- 
late the statististics of a magnetic system by generating a 
Markov chain of spin configurations and taking an average 
over the latter. Each step of this chain (a MC step) is a 
stochastic transition of the system from one state to an- 
other, subjected to the condition of the detailed balance. 
Usually a MC step consists in generating a new trial ori- 
entation of a spin vector on a lattice site i and calculating 
the ensuing energy change of the system AE. The trial 
confuguration is accepted as a new configuration if 



exp(-AE/T) > ft(0,l), 



(74) 



8 



H. Kachkachi and D. A.Garanin: Spin-wave theory for finite classical magnets and superparamagnetic relation 



m, M 



where 7?.(0, 1) is a random number in the interval (0, 1), 
otherwise the old configuration is kept. As follows from 
Eq. (|74|), for AE < the trial orientation is accepted with 
a probability 1. The trial orientation can be a completely 
random orientation, or a random orientation in the vicin- 
ity of the initial orientation of the spin , which is more 
appropriate at low temperatures. For the Ising model, the 
trial orientation is generated by a flip of Si with a probabil- 
ity 1/2. The MC steps are performed sequentially or ran- 
domly for all lattice sites. This conventional version of the 
MC method is not efficient for systems of finite size at low 
temperatures and small fields, if one is interested in the in- 
duced magnetization m. The Boltzmann distribution over 
the directions of the particle's magnetic moment M of Eq. 
(§) is achieved by rotations of M itself rather than by ro- 
tations of individual spins s^. Indeed, each spin s.; is acted 
upon by the strong exchange field H^j = J2j Jij s j ~ Jo, 
and in the typical case H -C Jo all trial configurations 
with the direction of s, significantly differing from that 
of its neighbors are rejected with a probability close to 1. 
Thus in the standard MC procedure directions of individ- 
ual spins can only change little by little, and the resulting 
change of M is extremely slow. For the Ising model the 
situation is even worse since the spin geometry is discrete 
and there are no small changes of spin directions, whereas 
a flip of a single spin against the exchange field has an ex- 
ponentially small probability. Hence if one starts in zero 
field with the configuration of all spins up or all spins 
down, the magnetization m will practically never relax to 
zero. This drawback can be remedied by augmenting the 
procedure by a global rotation (GL) of the particle's spins 
to a new trial direction of M and calculating the energy 
change. That is, before turning single spins on all lattice 
sites, one computes M, generates its new orientation M' 
and obtains the energy difference 



AE = A/"H • (M' - M). 



(75) 
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Fig. 2. Field dependence of the intrinsic magnetization M 
and the induced magnetization m of the Ising model on the sc 
lattice with fbc for different temperatures. 
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Fig. 3. Field dependence of the intrinsic magnetization M 
and the induced magnetization m of the Heisenberg model on 
the sc lattice with fbc for different temperatures. 



If the new orientation is accepted according to Eq. (|74|), 
one turns all spins by an appropriate angle and proceeds 
with the standard Metropolis method recapitulated above. 
In small fields (x = AfH/T < 1) relaxation of the induced 
magnetization m becomes much slower than that of the 
intrinsic magnetization M, and one needs much more MC 
steps to find the former than the latter with the same pre- 
cision. If in the procedure each global rotation is coupled 
with subsequent turning of single spins on all lattice sites 
i, making enough global rotations to achieve a required 
precision for m costs much more computer time for larger 
particle sizes. Thus a natural idea is to make many global 
rotations and gather the data for m after each GL before 
proceeding to the conventional (single-spin) part of the 
Metropolis algorithm. This improved method is especially 
fast for the isotropic Heisenberg or Ising models where 
the energy change is given by Eq. ( |75|) since, after M has 
been initially computed, each of its subsequent rotations 
and calculations of AE requires 0(1) operations. In con- 
trast, for systems with anisotropy one has to perform a 
sum over all lattice sites for each orientation of M, i.e., to 
make 0(Af) operations. 



Finally, we mention that for the Heisenberg model 5 3 
the running time of our programme with global rotations 
on a Pentium III/933 MHz is 160mn, for a precision of 
1CP 4 — 10~ 5 on the magnetisation. 

Figs. | and | show the results of our MC simulations 
for the Ising and Heisenberg models on a cubic lattice with 
size Af — 5 3 and free boundary conditions. The intrinsic 
magnetization M and induced magnetization m are plot- 
ted vs the scaled field x = AfH/T for different tempera- 
tures. We used the bulk Curie temperatures T c — 9 C T^ FA , 
where T c MFA — Jq/D is the mean- field Curie tempera- 
ture and 9 C is 0.751 for the Ising model and 0.722 for the 
Heisenberg model. One can see that the particle's mag- 
netic moment is aligned and thus m ~ M for x > 1, if 
T < T c . At T > T c the field alignes individual spins 
and this requires H > T, i.e., x > Af. The quadratic de- 
pendence of M(H) at small fields, which is phenomeno- 
logically described by Eq. (||) , manifests itself strongly at 
elevated temperatures. At low temperatures this depen- 
dence is much more difficult to see on the graph because 
the field-dependent part of M, which for the Heisenberg 
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Fig. 4. Scaled graph for the induced magnetization m of the 
Ising model on the sc lattice with fbc for different tempera- 
tures. Theoretical curves at low temperatures, B\(x) = tanhx, 
and at high temperatures, B ao (x), are shown by solid lines. 
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Fig. 5. Scaled graph for the induced magnetization m of 
the Heisenberg model on the sc lattice with fbc for different 
temperatures. Theoretical curves at low temperatures, Bs(x) — 
cothx — 1/a;, and at high temperatures, B 00 (x), are shown by 
solid lines. 



model is given by Eq. (^lj) , is for x ~ 1 proportional to a of 
Eq. (|54|), which is very small. For the Ising model there is 
practically no field dependence of M at low temperatures 
since M is very close to 1. The weak linear field depen- 
dence of M for the Heisenberg model which is visible on 
Fig. U at T = T c /4 will be quantitatively explained below. 

Figs.0 and [| show that the superparamagnetic relation 
of Eq. p) with the Langevin function Bd(x) in place of 
F(x) is a very good approximation everywhere below T c , 
for both Ising and Heisenberg models. On the other hand, 
above T c Eq. (|l|) with the function B oa (x) of Eq. is 
obeyed. The difference between these limiting expressions 
decreases with increasing the number D of spin compo- 
nents and disappears in the spherical limit (D — > oo). 

In Fig. | we compare theoretical predictions of Sec. 
| for the Heisenberg model at T = T c /4 with our MC 
results. For our small size M — 5 3 the square-root field 
dependence of the magnetization [the third line of Eq. 
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Fig. 6. Comparison of the theoretical and MC results for 
the field dependences of the magnetizations M and m for the 
Heisenberg model at T = T c /4. 



(pq)] does not arise and finite-size corrections are very 
important. For M one should use Eq. (|6l|), where t and 
a are given by Eq. (|54| ) with the numerically exact values 
W N = 1.99 and c N = 1.66 for the fbc model [cf. Eqs. © 
and @]. This yields t ~ 0.119 and a ~ 1.20 x 10" 4 . The 
corresponding theoretical dependence M(H) is practically 
a straight line which goes slightly above the MC points. 
This small discrepancy can be explained by the fact that 
the applicability criterion of our analytical method, t <S 1, 
is not strongly satisfied. For comparison we also plotted 
the theoretical M(H) for the unrealistic model with peri- 
odic boundary conditions. Here one has Wn = 1-25 and 
cat = 0.20, thus t ~ 0.075 and a ~ 1.45 x 10" 5 , so M(H) 
goes noticeably higher and with a much smaller slope. The 
quadratic field dependence of M in the region x < 1 is not 
seen at this low temperature since the value of a is very 
small and thus much more accurate MC simulations are 
needed. We have not performed these simulations because 
the corresponding effects are very small. We also plotted 
in Fig. [| the field dependence of m given by Eq. (^) in 
comparison with our MC data. The agreement is reason- 
ably good for m as well. 



5 Discussion 

We have performed analytical and numerical investigation 
of the magnetic field dependence of the intrinsic magneti- 
zation M and induced magnetization m of the Ising and 
isotropic classical Heisenberg models on cubic lattices of 
finite size. For the latter, we obtained explicit analytical 
results for both M(H,T) and m(H,T) at low tempera- 
tures with the help of a spin-wave theory singling out the 
global-rotation mode. These results are in accord with our 
MC simulation data. 

We investigated the validity of the superparamagnetic 
relation m(H, T) = M(H,T)B D [M{H,T)AfH/T], where 
Bjj(x) is the Langevin function and D is the number 
of spin components. Both general arguments of Sec. || 
and explicit low-temperature results for the Heisenberg 
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model show that this is not an exact relation for any fi- 
nite D. Nevertheless, it is an extremely good approxima- 
tion in the whole range below T c for not too small parti- 
cles, since, for the Heisenberg model, its error behaves as 
[77 ( J N)] 2 , where N is the linear particle size. For N 3> 1, 
a crossover to the high-temperature form of the relation 
above, which utilizes the Langevin function of the spher- 
ical model Bqo^x), occurs in a close vicinity of T c . The 
difference between the low- and high-temperature forms 
of the superparamagnetic relation decreases with D and 
disappears in the spherical limit, rendering this relation 
exact 

D. A. Garanin is endebted to the CNRS and Labo- 
ratoire de Magnetisme et d'Optique for the warm hos- 
pitability extended to him during his stay in Versailles in 
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